clear; close all; clc;


want_durs = [10 0.15];
freqs= [300 600 1500 3000]; %just a random frequency

subjects = {'003','008','009','010','011','012','015'}; % we chose several test subjects.
% subjects = {'003','008'}; % we chose several test subjects.

plot_each =0;
projection1 = 0; % neariest center
% projection1 = 1; % first project to the unit circle
Num_centers = 1;

true_angle = [-80 -65 -55 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 55 65 80];
Elevation_number = 9;

f1=figure('position',[100 70 1100 520]);
pause(0.1)
sig='* ';
figss = reshape([ 1 2 3 4;5 6 7 8],1,8)
TAs = 'AEBFCGDH';
c_fig = 0;
for ifreql = 1:length(freqs)
    freq = freqs(ifreql)
    for idur = 1:length(want_durs)
        c_fig=c_fig+1;
        [y,fs]=audioread('F_us_8k.wav');
        y = y(fs:end,1); % Always truncate out the 1st 1 second because it has no sound.
        want_dur=want_durs(idur)
        y = y(1:round(want_dur*fs),1);
        ramplength = 5e-3;
        wholeramp=f_getRamp(ramplength,want_dur,fs);
        if length(y)>length(wholeramp)
            wholeramp=[wholeramp zeros(1,length(y)-length(wholeramp))];
        end
        y = y.*wholeramp';
        %         sound_duration = length(y)/fs
        %--------------------------------------------------------------
        % RIR info
        cline = 'rbmkgcy';
        for isub = 1:length(subjects)
            eval(['load hrir_final' subjects{isub} '.mat']);
            for iangle = 1:length(true_angle)
                H2l=squeeze(hrir_l(iangle,Elevation_number,:));  H2r=squeeze(hrir_r(iangle,Elevation_number,:));
                s2l=filter(H2l,1,y);
                s2r=filter(H2r,1,y);
                x1=s2l;  %left channel
                x2=s2r;
                % short time fourier transform:
                [SL,SR,freq_pos,flag]=mySTFT(x1,x2,fs);
                % normalization:
                [SL1,SR1]=Normalize(SL,SR);
                [Y Ind] = min(abs(freq_pos - freq));
                X_P = [real(SR1(Ind,flag>0))' imag(SR1(Ind,flag>0))'];
                deg_sign=char(0176);
                [IDX, C] = kmeans(X_P, Num_centers,'Distance','cityblock','Replicates',5);
                if plot_each == 1
                    openfig('model2.fig'); hold on
                    C
                    plot(C(1,1),C(1,2),'g*','Markersize',24,'LineWidth',3);hold on;
                    if Num_centers > 1
                        plot(C(2,1),C(2,2),'ro','Markersize',24,'LineWidth',3);hold on;
                    end
                    if Num_centers > 2
                        plot(C(3,1),C(3,2),'mo','Markersize',24,'LineWidth',3);hold on;
                    end
                    grid on
                    set(gca,'Xtick',-1:0.5:1)
                    set(gca,'Ytick',-1:0.5:1)
                    axis([-1 1 -1 1])
                    hold on;
                    plot(-1:1,[0,0,0],'k', 'LineWidth', 0.5);
                    plot([0,0,0],-1:1,'k','LineWidth',0.5)
                    set(gca,'fontsize',12)
                    title([tps{roomtype} ', ' num2str(true_angle(iangle)) deg_sign ],'fontsize',15,'fontweight','bold')
                    xlabel('Real','fontsize',15)
                    ylabel('Imaginary','fontsize',15)
                    hold on
                    % now plot the line:
                    C_new = C./sqrt(C(1)^2+C(2)^2);
                    plot([0 C_new(1)],[0 C_new(2)],'r');
                    axis square
                end
                clear x1 y1 y0 m
                eval(['load Freq_centers' num2str(freq) ])
                out = f_angle_class_center(C,F_centers,projection1); % cluster centers
                classified_angle(isub,iangle) = true_angle(out.class);
                error(isub,iangle) = abs(true_angle(iangle)-classified_angle(isub,iangle));

            end % iangle
            figure(f1)
            subplot(2,4,figss(c_fig))
            s=plot(true_angle,classified_angle(isub,:),'*-');set(s,'color',cline(isub));hold on;
            axis square
            xlabel('True Angle (^{\circ})','fontsize',9);ylabel('Classified Angle (^{\circ})','fontsize',9);
            set(gca,'fontsize',8,'xtick',-100:50:100)

            pause(0.1)

        end
        % r_corr
        all_1 = 0.1.*round(10*mean(mean(error)));
        text(-95,95,['Mean {\it e}=' num2str(all_1) '^{\circ}']);
        for i=1:length(subjects)
            legendsall{i}=['S' subjects{i} ', {\it e}=' num2str(mean(error(i,:))) '^{\circ}'];
        end
        hold on;plot(-90:90,-90:90,'k:');
        box off;
        s=title([TAs(c_fig) '. ' num2str(freq) ' Hz, ' num2str(want_dur) ' s']);set(s,'fontsize',12)
    end % dur

    [h,p]=ttest(error(1,:),error(2,:));
    pause(0.1)
    p_all(ifreql)=p;
end % freq
p_all
exportgraphics(gcf, 'Figure_9_duration.jpg', 'Resolution', 300,'ContentType', 'vector')